set more off

* the data used in this paper may not be publicly released and aspects of the data cleaning files may violate terms of our data use agreeements.
* our code takes a final analysis dataset as given, and provides code of all actual analyses in the paper.
* this file reflects analysis for period 2 (2006-2012), bottled water, batteries, flashlights ... for HOUSEHOLD-LEVEL DATA.


*************************************************************************************************************************************
* ANALYSIS CONSTRUCTION - HH level data
*************************************************************************************************************************************

use period3_county_treatments, clear
tsset, clear

merge 1:m countyid statadate using period3_panel_expenditures_grocery
* merge 1:1 countyid statadate using period3_county_expenditures_grocery_temp
	* merge=1 is counties without expenditures in our dataset
	* merge=2 is month without full info ... jan 2004, dec 2012 ... or one of the dropped counties
keep if _merge==3
drop _merge

destring countyid, gen(countyid_num)

tsset household_code statadate


*********************************************************
* event study

* note that these are point-wise confidence intervals

svyset [pw=projection_factor],psu(household_code)

* event study all goods unconditional

foreach z in "0000" {

preserve

svy:mean expend_`z'
gen mean=_b[expend_`z']
replace expend_`z'=expend_`z'-mean
* note since indiv data do not rescale by dividing by 1000 as in store-level event study graphs

local win = 15 
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}

*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_* projection_factor household_code 

forvalues x=0/`win' {
	svy:mean expend_`z'_l`x'
	gen expend3_`z'_l`x'=_b[expend_`z'_l`x']
	gen expend2_`z'_l`x'=_se[expend_`z'_l`x']
	replace expend_`z'_l`x'=expend3_`z'_l`x'
	drop expend3*
}

forvalues x=1/`win' {
	svy:mean expend_`z'_f`x'
	gen expend3_`z'_f`x'=_b[expend_`z'_f`x']
	gen expend2_`z'_f`x'=_se[expend_`z'_f`x']
	replace expend_`z'_f`x'=expend3_`z'_f`x'
	drop expend3*
}

collapse (mean) expend_`z'_* expend2_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi -4 0 2.5 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(-4(2)2) xlabel(-15(3)15) xtitle("Days Before and After Storm Landfall") legend(off) title("Mean expenditures - all tracked goods") saving(Fig4A_w, replace)

restore

}



* event study emergency supply goods conditional
foreach z in "1487" {

preserve

replace expend_`z'=. if expend_0000==0

svy:mean expend_`z'
gen mean=_b[expend_`z']
replace expend_`z'=expend_`z'-mean
* note since indiv data do not rescale by dividing by 1000 as in store-level event study graphs

local win = 15 
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}


*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_* projection_factor household_code 

forvalues x=0/`win' {
	svy:mean expend_`z'_l`x'
	gen expend3_`z'_l`x'=_b[expend_`z'_l`x']
	gen expend2_`z'_l`x'=_se[expend_`z'_l`x']
	replace expend_`z'_l`x'=expend3_`z'_l`x'
	drop expend3*
}


forvalues x=1/`win' {
	svy:mean expend_`z'_f`x'
	gen expend3_`z'_f`x'=_b[expend_`z'_f`x']
	gen expend2_`z'_f`x'=_se[expend_`z'_f`x']
	replace expend_`z'_f`x'=expend3_`z'_f`x'
	drop expend3*
}

collapse (mean) expend_`z'_* expend2_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi -0.2 0 1.3 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(0(.4)1.2) xlabel(-15(3)15) xtitle("Days Before and After Storm Landfall") legend(off) title("Per household expenditures on bottled water, 2004-2012") saving(Fig3A_w, replace)

restore

}



* event study emergency supply goods conditional
foreach z in "7870" {

preserve

replace expend_`z'=. if expend_0000==0

svy:mean expend_`z'
gen mean=_b[expend_`z']
replace expend_`z'=expend_`z'-mean
* note since indiv data do not rescale by dividing by 1000 as in store-level event study graphs

local win = 15 
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}


*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_* projection_factor household_code 

forvalues x=0/`win' {
	svy:mean expend_`z'_l`x'
	gen expend3_`z'_l`x'=_b[expend_`z'_l`x']
	gen expend2_`z'_l`x'=_se[expend_`z'_l`x']
	replace expend_`z'_l`x'=expend3_`z'_l`x'
	drop expend3*
}


forvalues x=1/`win' {
	svy:mean expend_`z'_f`x'
	gen expend3_`z'_f`x'=_b[expend_`z'_f`x']
	gen expend2_`z'_f`x'=_se[expend_`z'_f`x']
	replace expend_`z'_f`x'=expend3_`z'_f`x'
	drop expend3*
}

collapse (mean) expend_`z'_* expend2_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi -0.1 0 1.0 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(0(.3).9) xlabel(-15(3)15) xtitle("Days Before and After Storm Landfall") legend(off) title("Per household expenditures on batteries, 2004-2012") saving(Fig3B_w, replace)

restore

}


* event study emergency supply goods conditional
foreach z in "7871" {

preserve

replace expend_`z'=. if expend_0000==0

svy:mean expend_`z'
gen mean=_b[expend_`z']
replace expend_`z'=expend_`z'-mean
* note since indiv data do not rescale by dividing by 1000 as in store-level event study graphs

local win = 15 
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}


*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_* projection_factor household_code 

forvalues x=0/`win' {
	svy:mean expend_`z'_l`x'
	gen expend3_`z'_l`x'=_b[expend_`z'_l`x']
	gen expend2_`z'_l`x'=_se[expend_`z'_l`x']
	replace expend_`z'_l`x'=expend3_`z'_l`x'
	drop expend3*
}


forvalues x=1/`win' {
	svy:mean expend_`z'_f`x'
	gen expend3_`z'_f`x'=_b[expend_`z'_f`x']
	gen expend2_`z'_f`x'=_se[expend_`z'_f`x']
	replace expend_`z'_f`x'=expend3_`z'_f`x'
	drop expend3*
}

collapse (mean) expend_`z'_* expend2_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi -0.1 0 0.33 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(0(.1)0.3) xlabel(-15(3)15) xtitle("Days Before and After Storm Landfall") legend(off) title("Per household expenditures on flashlights, 2004-2012") saving(Fig3C_w, replace)

restore

}




* event study - no trips

preserve

sum fillin
gen mean=r(mean)
replace fillin=fillin-mean

local win = 15 
forvalues x=0/`win' {
	gen fillin_l`x'=l`x'.fillin
}
forvalues x=1/`win' {
	gen fillin_f`x'=f`x'.fillin
}

*keep if hlandfall_false==1
keep if hlandfall==1

keep fillin_* projection_factor household_code
 
forvalues x=0/`win' {
	svy:mean fillin_l`x'
	gen fillin3_l`x'=_b[fillin_l`x']
	gen fillin2_l`x'=_se[fillin_l`x']
	replace fillin_l`x'=fillin3_l`x'
	drop fillin3*
}


forvalues x=1/`win' {
	svy:mean fillin_f`x'
	gen fillin3_f`x'=_b[fillin_f`x']
	gen fillin2_f`x'=_se[fillin_f`x']
	replace fillin_f`x'=fillin3_f`x'
	drop fillin3*
}


collapse (mean) fillin_* fillin2_* 

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename fillin_f`x' fillin_`y'
	rename fillin2_f`x' fillin2_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename fillin_l`x' fillin_`y'
	rename fillin2_l`x' fillin2_`y'
}

gen id=1
reshape long fillin_ fillin2_, i(id) j(time)
drop id
replace time=time-50

gen ci=fillin2_*1.96
gen ul=fillin_+ci
gen ll=fillin_-ci

twoway (rarea ul ll time, color(gs14)) (connected fillin_ time, msize(small)) (pcarrowi -0.05 0 0.1 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(-0.05(.05)0.1) xlabel(-15(3)15) xtitle("Days Before and After Storm Landfall") legend(off) title("% households recording no trips, 2004-2012") saving(Fig4B_w, replace)

restore



****************************************************************************************************
* ADDITIONAL PROCESSING FOR REGRESSION ANALYSIS
****************************************************************************************************

*merge in weather. note some winter days are missing weather data. we do not use winter days, so not an issue. 
merge m:1 countyid statadate using county_weather_nearest
keep if _merge==3
drop _merge

* create year dummies
tab year, gen(years)

* keep only months in hurricane seasons and create month dummies
gen month=month(statadate)
keep if inrange(month,6,11)
tab month, gen(months)

* generate logged depvars - NOTE ADD $1 and log
foreach code in "0000" "1487" "7870" "7871" {
	gen lexpend_`code'=log(expend_`code'+1)
}

* drop 335 HHs with no tracked grocery expenditures during our period.
bys household_code: egen temp=max(expend_0000)
drop if temp==0
drop temp


***********************************************************************************************************************************
* REGRESSIONS
***********************************************************************************************************************************

* main regressions - with HH level fixed effects

log using period3_results_HHFEs, replace text

* unconditional

foreach code in "0000" "1487" "7870" "7871" {
	areg lexpend_`code' htreat_before hlandfall htreat_after prcp tmax tmin months1-months5 years1-years8 [pweight=projection_factor], a(household_code) vce(cluster countyid_num)
}

foreach code in "0000" "1487" "7870" "7871" {
	areg lexpend_`code' htreat_before96 htreat_before72 htreat_before48 htreat_before24 hlandfall htreat_after prcp tmax tmin months1-months5 years1-years8 [pweight=projection_factor], a(household_code) vce(cluster countyid_num)
}

* percent no trips

areg fillin htreat_before hlandfall htreat_after prcp tmax tmin months1-months5 years1-years8 [pweight=projection_factor], a(countyid_num) vce(cluster countyid_num)

areg fillin htreat_before96 htreat_before72 htreat_before48 htreat_before24 hlandfall htreat_after prcp tmax tmin months1-months5 years1-years8 [pweight=projection_factor], a(household_code) vce(cluster countyid_num)

* conditional on consuming

foreach code in "0000" "1487" "7870" "7871" {
	areg lexpend_`code' htreat_before hlandfall htreat_after prcp tmax tmin months1-months5 years1-years8 if expend_0000!=0 [pweight=projection_factor], a(household_code) vce(cluster countyid_num)
}

foreach code in "0000" "1487" "7870" "7871" {
	areg lexpend_`code' htreat_before96 htreat_before72 htreat_before48 htreat_before24 hlandfall htreat_after prcp tmax tmin months1-months5 years1-years8 if expend_0000!=0 [pweight=projection_factor], a(household_code) vce(cluster countyid_num)
}

STOP
STOP
STOP
STOP

